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Abstract. In this paper, we investigate interpolatory projection framework for model reduction 
of descriptor systems. With a simple numerical example, we first illustrate that employing sub- 
space conditions from the standard state space settings to descriptor systems generically leads to 
unbounded H2 or Woo errors due to the mismatch of the polynomial parts of the full and reduced- 
order transfer functions. We then develop modified interpolatory subspace conditions based on the 
deflating subspaces that guarantee a bounded error. For the special cases of index-1 and index-2 de- 
scriptor systems, we also show how to avoid computing these deflating subspaces explicitly while still 
enforcing interpolation. The question of how to choose interpolation points optimally naturally arises 
as in the standard state space setting. We answer this question in the framework of the %2-norm by 
extending the Iterative Rational Krylov Algorithm (IRKA) to descriptor systems. Several numerical 
examples are used to illustrate the theoretical discussion. 
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1. Introduction. We discuss interpolatory model reduction of differential- 
algebraic equations (DAEs), or descriptor systems, given by 

Ex(i) = Ax(t)+Bu(t), . . 

y(t) = Cx(tj+Du(i), [LA) 

where x(t) G E™, u(i) G R m and y(t) G R p are the states, inputs and outputs, 
respectively, E e R" x " is a singular matrix, A G R nxn , B G R nxm , C G R px, \ 



and D G K pxm . Taking the Laplace transformation of system (1.1) with zero initial 
condition x(0) = 0, we obtain y(s) = G(s)u(s), where u(s) and y(s) denote the 
Laplace transforms of u(t) and y(i), respectively, and G(s) = C(sE — A) _1 B + D 
is a transfer function of (JOJ. By following the standard abuse of notation, we will 
denote both the dynamical system and its transfer function by G. 



Systems of the form (1.1) with extremely large state space dimension n arise in 
various applications such as electrical circuit simulations, multibody dynamics, or 
semidiscretized partial differential equations. Simulation and control in these large- 
scale settings is a huge computational burden. Efficient model utilization becomes 
crucial where model reduction offers a remedy. The goal of model reduction is to 



replace the original dynamics in (1.1) by a model of the same form but with much 
smaller state space dimension such that this reduced model is a high fidelity approx- 
imation to the original one. Hence, we seek a reduced-order model 

Ex(t) = Ax(i)+Bu(i), 
y(f) = Cx(t)+Du(t), 1 • ' 

where E, A G R rxr , B G R rxm , C G R pxr , and D G R pxm such that r < n, and the 
error y — y is small with respect to a specific norm over a wide range of inputs u(t) 
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wit h bo unded energy. In the frequency domain, this means that the transfer function 



of (1.2 1 given by G(s) = C(sE — A) 1 B + D approximates G(s) well, i.e., the error 



G(s) — G(s) is small in a certain system norm. 



The reduced-order model (1.2 1 can be obtained via projection as follows. We 
first construct two n x r matrices V and W, approximate the full-order state x(i) by 
Vx(i), and then enforce the Petrov-Galerkin condition 

W T f EVx(i) - AVx(i) - B u(t)) = 0, y(<) = CVx(i) + Du(t). 



As a result, we obtain the reduced-order model (1.2 1 with the system matrices 



E = W T EV, A = W T AV, 

(1.3) 

B = W B, C = CV, D = D. 

The projection matrices V and W determine the subspaces of interest and can be 
computed in many different ways. 

In this paper, we consider projection-based interpolatory model reduction me- 
thods, where the choice of V and W enforces certain tangential interpolation of the 
original transfer function. These methods will be presented in Section [2] in more de- 
tail. Projection-based interpolation with multiple interpolation points was initially 
proposed by Skelton et. al. in |3TJ [35] ■ Grimme [TU] has later developed a numeri- 
cally efficient framework using the rational Krylov subspace method of Ruhe [21] . The 
tangential rational interpolation framework, we will be using here, is due to a recent 
work by Gallivan et al. [9]. 

Unfortunately, it is often assumed that extending interpolatory model reduction 
from standard state space systems with E = I to descriptor systems with singular 
E is as simple as replacing I by E. In Section [2] we present an example showing 
that this naive approach may lead to a poor approximation with an unbounded error 
G(s) — G(s) although the classical interpolatory subspace conditions are satisfied. 
In Section |3j we modify these conditions in order to enforce bounded error. The 
theoretical result will take advantage of the spectral projectors. Then using the new 
subspace conditions, we extend in Section [4] the optimal %2 model reduction method 
of [TS] to descriptor systems. Sections [3] and [4] make explicit usage of deflating sub- 
spaces which could be numerically demanding for general problems. Thus, for the 
special cases of index-1 and index-2 descriptor systems, we show in Sections [5] and 
[6] respectively, how to apply interpolatory model reduction without explicitly com- 
puting the deflating subspaces. Theoretical discussion will be supported by several 
numerical examples. In particular, in Section |5.2[ we present an example, where the 
balanced truncation approach |26j is prone to failing due to problems solving the gen- 
eralized Lyapunov equations, while the (optimal) interpolatory model reduction can 
be effectively applied. 

2. Model reduction by tangential rational interpolation. The goal of 
mod el reduction by tangential interpolation is to construct a reduced-order model 



(1.2) such that its transfer function G(s) interpolates the original one, G(s), at se- 
lected points in the complex plane along selected directions. We will use the notation 
of pQ to define this problem more precisely: Given G(s) = C(sE — A) _1 B + D, the 
left interpolation points {Hi} q i=1 , Hi € C, together with the left tangential directions 
{ci}* =1 , Ci £ C p , and the right interpolation points {<7j}^ =1 , <Jj £ C, together with the 
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right tangential directions {bj}j =1 , bj E C m , we seek to find a reduced-order model 
G(s) = C(sE — A) _1 B + D that is a tangential interpolant to G(s), i.e., 



cfG( Mi )=cfG( M i) ! 



3 



1,. 
1, 



,q, 



(2.1) 



Through out the paper, we will assume q — r, meaning that the same number of left 
and right interpolation points are used. In addition to interpolating G(s), one might 
ask for matching the higher-order derivatives of G(s) along the tangential directions 
as well. This scenario will also be handled. 

By combining the projection-based reduced-order modeling technique with the 
interpolation framework, we want to find the n x r matrices W and V such that 



the reduced-order model (1.2 1, (1.3| satisfies the tangential interpolation conditions 



(2.1 ). This approach is called projection-based interpolatory model reduction. How to 



enforce the interpolation conditions via projection is shown in the following theorem, 
where the ^-th derivative of G(s) with respect to s evaluated at s — a is denoted by 
GW( ff ). 

Theorem 2.1. [TJ [5] Let a, \i e C be such that sE — A and sE — A are both 
invertible for s = ex, fi, and let b € C m and c € C p be fixed nontrivial vectors. 
1. If 



((crE- A) _1 E) J 1 (crE- A) _1 Bb e Im(V), j = 1, 

then G< £ )(cr)b = G^(a)b for £ = 0, 1, . . . , N - 1. 
2. If 



(2.2) 



((/iE- Ar T E T ) J 1 (^E-A)" r C T ceIm(W), j = l,...,M, (2.3) 

then c T G' ( V) =c t GM(;i) for I = 0, 1, . . . , M - 1. 
3. If both Q and Q hold, and if a = fx, then c T G^{a)b = c T G^(a)b for 
e = 0,l,...,M + N + l. 
One can see that to solve the rational tangential interpolation problem via pro- 
jection all one has to do is to construct the matrices V and W as in Theorem |2.1| 
The dominant cost is to solve sparse linear systems. We also note that in Theorem |2.1| 
the values that are interpolated are never explicitly computed. This is crucial since 
that computation is known to be poorly conditioned [8]. 



To illustrate the result of Theorem 2.1 for a special case of Hcrmitc bi-tangential 
interpolation, we take the same right and left interpolation points {ci}J" =1 , left tangen- 
tial directions {q}[ =1 , and right tangential directions {bi}[ =1 . Then for the projection 
matrices 



V= [(cr 1 E-A)- 1 Bb 1 , (oyE — A) _1 Bb r ] , 

W = [(o- x E - A)- T C T ci, • • • , (cr r E - A)- T C T c r ] 



the reduced-order model G(s) = C(sE — A) 1 B + D as in ( 1.3 1 satisfies 
Gfajbi = Gfajbi, cfGfa) = cfG(iTj), cj G'^)^ = cf G'^)^ 



(2.4) 
(2.5) 



(2.6) 



for i = 1, • • • , r, provided cr^E — A and <t,-E — A are both nonsingular. 
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Note that Theorem |2.1| does not distinguish between the singular E case and the 
standard state space case with E = I. In other words, the interpolation conditions 
hold regardless as long as the matrices <7;E — A and <t,E — A are invertible. This is 
the precise reason why it is often assumed that extending interpolatory-based model 
reduction from G(s) = C(sl- A) _1 B + D to G(s) = C(sE- A) _1 B + D is as simple 
as replacing I by E. However, as the following example shows, this is not the case. 

Example 2.1. Consider an RLC circuit modeled by an index-2 SISO descriptor 
system (1.1 1 of order n — 765 (see, e.g., 20J for a definition of index). We approxi- 
mate this system with a model ( 1.2 ) of order r — 20 using Hermite interpolation. The 



carefully chosen interpolation points were taken as the mirror images of the dominant 
poles of G(s). Since these interpolation points are known to be good points for model 
reduction [11[ 113] , one would expect the interpolant to be a good approximation as 
well. However, the situation is indeed the opposite. Figure |2.1| shows the amplitude 
plots of the frequency responses G(iuj) and G(uS) (upper plot) and that of the er- 
ror G(zw) — G(iui) (lower plot). One can see that the error G(iu>) — G(ilu) grows 
unbounded as the frequency ui increases, and, hence, the approximation is extremely 
poor with unbounded H2 and Hoo error norms even though it satisfies Hermite inter- 
polation at carefully selected effective interpolation points. 



Amplitude Bode plots of G and G 




Amplitude Bode plot of the error G — G 




Fig. 2.1. Example 2.1' amplitude plots o/G(zw) and G(tw) (upper); the absolute error 
\G(iw) - G(jw)| (lower). 



The reason is simple. Even though E is singular, E = W T EV will generically 
be a nonsingular matrix assumin g r < rank(E). In this case, the transfer function 
G(s) of the reduced-order model (1.2 1 is proper, i.e., lim G(s) < 00, although G(s) 



might be improper. Hence, the special care needs to be taken in order to match the 
polynomial part of G(s). We note that the polynomial part of G(s) has to match 
that of G(s) exactly. Otherwise, regardless of how good the interpolation points are, 
the error will always grow unbounded. For the very special descriptor systems with 
the proper transfer functions and only for interpolation around s — 00, a solution is 
offered in [5]. For descriptor systems of index 1, where the polynomial part of G(s) 
is a constant matrix, a remedy is also suggested in [1] by an appropriate choice of D. 
However, the general case is remained unsolved. We will tackle precisely this problem, 
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where is a descriptor system of higher index, its transfer function G(s) may have 
a higher order polynomial part and interpolation is at arbitrary points in the complex 
plane. Thereby, the spectral projectors onto the left and right deflating subspaces of 
the pencil AE — A corresponding to finite eigenvalues will play a vital role. Moreover, 
we will show how to choose interpolation points and tangential directions optimally 
for interpolatory model reduction of descriptor systems. 

3. Interpolatory projection methods for descriptor systems. As stated 
above, in order to have bounded "Hoo and T~Li errors, the polynomial part of G(s) has 
to match the polynomial part of G(s) exactly. Let G(s) be additively decomposed as 

G( s ) = G sp ( s )+P( s ), (3.1) 

where G sp (s) and P(s) denote, respectively, the strictly proper part and the polyno- 
mial part of G(s). We enforce the reduced-order model G(s) to have the decomposi- 
tion 

G( s ) = G sp (s) + P( s ) (3.2) 

with P(s) = P(s). This implies that the error transfer function does not contain 
a polynomial part, i.e., 

G err (s) = G(s) - G(s) = G sp (s) - G sp (s) 

is strictly proper meaning lim G err (s) = 0. Hence, by making G sp (s) to interpolate 

G sp (s), we will be able to enforce that G(s) interpolates G(s). This will lead to 
the following construction of G(s). Given G(s), we create W and V satisfying new 
subspace conditions such that the reduced-order model G(s) obtained by projection as 



in ( 1.3 ) will not only satisfy the interpolation conditions but also match the polynomial 
part of G(s). 

Theorem 3.1. Given a full-order model G(s) — C(sE — A) _1 B + D, define 
P; and P r to be the spectral projectors onto the left and right deflating subspaces of 
the pencil AE — A corresponding to the finite eigenvalues. Let the columns of 
and Vqo span the left and right deflating subspaces of AE — A corresponding to the 
eigenvalue at infinity. Let a . fi € C be interpolation points such that sE — A and 
sE — A are nonsingular for s = a,fi, and let b € C m and c£p. Define *Vf and W/ 
such that 

Im(V/) =span|((crE- A)" 1 E)^ 1 (crE-A)" 1 P ; Bb, j = l,...,ivj, (3.3) 
Im(W» =span{((^E-A)- T E T ) : ' _1 ( ) uE-A)- T P^C T c, j = 1, m} . (3.4) 
Then with the choice of W = [W/, Woo] and V = [V/, Vqo], the reduced-order 



model G(s) = C(sE — A) 1 B + D obtained via projection as in (1.3) satisfies 

1. P(s) = P(s) L 

2. G^(a)b = G^(cr)b for i = 0, 1, . . . , N - 1, 

3. c T GW(/i) = c T G^(Ai) for I = 0, 1, „ . , M - 1. 

If a = ii, we have, additionally, c T G^\a)b = c T G^(cr)b for £ = 0, M + N + l. 
Proof. Let the pencil AE — A be transformed into the Weierstrass canonical form 



E = S 



In, 
N 



T"\ A = S 



J 

L_ 



T , (3.5) 
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where S and T are nonsingular and N is nilpotcnt. Then the projectors P/ and P T 
can be represented as 



















(3.6) 



Let T = [Ti, T 2 ] and S 1 = [Si, S 2 ] T be partitioned according to E and A 
(3.5). Then the matrices and Voo take the form 

(I-PfjWoo, Voo =T 2 R T = (I-P r )V 00 



w.. 



S 2 R 



2*VS 



with nonsingular R5 and Ry. Furthermore, the strictly proper and polynomial parts 
of G(s) in (3.1) are given by 

G sp (s) = CT^I^- JJ-^fB, 
and P(s) = CT 2 (sN-I„ Do )- 1 S 2 n B + D, 



respectively. It follows from (13.51) and (3.6) that 



EP r = P;E, 

(sE- A)- 1 P ( 



AP r - P,A, 

P r (aE-A)- x , 



and, hence, 

W f = PfW f and V/ = P r V> . 
Then the system matrices of the reduced-order model have the form 



(3.7) 



E = W T EV 



A = W AV = 



B = W J B 



W/EV, 



WJEV 



WlEV f W^EV^ 

" VV/ AV, W/ AV X 

WlAV f W^AV^ 
W/B 



W/EV, 



WJA V/ 





W T EV 


W T AV 



C = CV = [CV), CVoo], 



D = D. 



W^B 

Thus, the strictly proper and polynomial parts of G(s) are given by 

G sp (s) - CV / (sW / T EV / - W/AV / )- 1 W/B, 
P(s) = CV 00 ( S W^EV 00 -W^AV CO )- 1 W^B + D 
= CT 2 (.sI - J)- 1 SjB + D = P(s). 

One can see that the polynomial parts of G(s) and G(s) coincide, and the proof of 
the interpolation result reduces to proving the interpolation cond ition s for the s trictly 
proper parts of G(s) and G(s). To prove this, we first note that (3.5 ) and (3.6) imply 
that 



CP r (crE — A) _1 P;B 



I 





crl J 




CT 

u u 

CT 1 ((7l-J)- 1 SrB = G sp ( ( r). 






-1 


I 


" 


o-N-I 
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Furthermore, it follows from the relations (3.7| that 



CP r V/ 



cv 



WjP z B = W/B. 



Due to the definitions of Vy and Wj in (3.3 1 and (3.4 1, respectively, Theorem 
gives 



2.1 



G 



sp (a)b = CPrVfiaWjEVj - W/AV / )" 1 W / r P;Bb = G sp (a)b, 



c J G sp ( M ) = c 1 CP r V f (jiWj EV f - Wf AV^Wf P,B = c 1 G 8p ( M ). 



Since both parts 1 and 2 of Theorem 2.1 hold, we have c T Gg p (<r)b = c T G' ap (a)b for 
a - = fj,. The other interpolatory relations for the derivatives of the transfer function 
can be proved analogously. □ 



Next, we illustrate that even though Theorem |3.1| has a very similar structure to 
that of Theorem |2.1 1 the saddle difference between these two results makes a big diffe- 
rence in the resulting reduced-order model. Towards this goal, we revisit Example |2.1| 
We reduce the same full-order model using the same interpolation points, but imposing 



the subspace conditions of Theorem |3.1[ instead. Figure 3T depicts the resulting 
amplitude plots of G(uS) and G(iw) (upper plot) and that of the error G(uj) — G(zw) 



(lower plot) whe n th e new subspace conditions of Theorem 3.1 are used. Unlike the 
case in Example 2.1 where the error G(tui) — G(ilu) grew unbounded, for the new 
reduced-order model, the maximum error is below 10 -2 and the error decays to zero 
as w approaches oo, since the polynomial part is captured exactly. 

Amplitude Bode plots of G and G 




10 



10 10 10 10 

Amplitude Bode plot of G - G 




10 



10 10 

freq (rad/sec) 



10 



10 



Fig. 3.1. Amplitude plots ofG(uS) andG(iuj) (upper); the absolute error \G{iui)~G(iuj) 
(lower). 



In some applications, the deflating subspaces of AE — A corresponding to the 
eigenvalues at infinity may have large dimension n^. However, the order of the system 
can still be reduced if it contains states that are uncontrollable and unobservable at 
infinity. Such states can be removed from the system without changing its transfer 
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function and, hence, preserving the interpolation conditions as in Theorem |3.1| In 
this case the projection matrices and can be determined as proposed in [25] 
by solving the projected discrete-time Lyapunov equations 

AXA T - EXE T = -(I - P/)BB T (I - P ; ) T , X = (I - P r )X(I - P r ) T , (3.8) 
A T YA - E T YE = -(I - P r ) T C T C(I - P r ), Y = (I — P;) T Y(I - P/). (3.9) 

Let Xc and Yc be the Cholesky factors of X — X c Xg and Y = Y c Yp, respectively, 
and let Y^AX C = [Ui, U ]diag(S, 0)[Vi, V ] T be singular value decomposition, 
where [Ui, Uo] and [Vi, Vo] are orthogonal and £ is nonsingular. Then the projec- 
tion matrices and can be taken as = YpUi and Voo = XcVi. Note 
that the Cholesky factors Xc and Yc can be computed directly using the generalized 
Smith method [27] . In this method, it is required to solve v — 1 linear systems only, 
where v is the index of the pencil AE — A or, equivalently, the nilpotence index of N 



in (3.5 1. The computation of the projectors P; and P r is, in general, a difficult prob- 
lem. However, for some structured problems arising in circuit simulation, multibody 
systems and computational fluid dynamics, these projectors can be constructed in 
explicit form that significantly reduces the computational complexity of the method; 
see [27] for details. 

4. Interpolatory optimal H2 model reduction for descriptor systems. 

The choice of interpolation points and tangential directions is the central issue in 
interpolatory model reduction. This choice determines whether the reduced-order 
model is high fidelity or not. Until recently, selection of interpolation points was 
largely ad hoc and required several model reduction attempts to arrive at a reasonable 
approximation. However, Gugercin et al. [15 introduced an interpolatory model 
reduction method for generating a reduced model G of order r which is an optimal 
%2 approximation to the original system G in the sense that it minimizes T^-norm 
error, i.e., 

||G-G|| Wa = mm ||G-G,.||„ 2 , (4.1) 

dim(G r )=r 



where 



I /.+« \ 1/2 

-L / ,1 „ , ^ . 1 2 



I g IIh 2 : =(^/ IIgMIIf^) (4.2) 



— 00 



and || ■ || p denotes the Frobenius matrix norm. Since this is a non-convex optimization 
problem, the computation of a global minimizer is a very difficult task. Hence, instead, 
one tries to find high-fidelity reduced models that satisfy first-order necessary opti- 
mally conditions. There exist, in general, two approaches for solving this problem. 
These are Lyapunov-based optimal H2 methods presented in [m|18j[25j ^> ^3] an d 
interpolation-based optimal H2 methods considered in [3] [4] |6j H2] HU [15j HH [23l [28] . 
While the Lyapunov-based approaches require solving a series of Lyapunov equations, 
which becomes costly and sometimes intractable in large-scale settings, the interpola- 
tory approaches only require solving a series of sparse linear systems and have proved 
to be numerically very effective. Moreover, as shown in [15], both frameworks are the- 
oretically equivalent that further motivates the usage of interpolatory model reduction 
techniques for the optimal H.2 approximation. 

For SISO systems, interpolation-based H2 optimality conditions were originally 
developed by Meier and Luenberger [55] . Then, based on these conditions, an effective 
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algorithm for intcrpolatory optimal Ti.2 approximation, called the Iterative Rational 
Krylov Algorithm (IRKA), was introduced in [T2JQ3]. This algorithm has also been 
recently extended to MIMO systems using the tangential interpolation framework, 
see E3 HE] for more details. 

The model reduction methods mentioned above, however, only deals with the 
system G(s) = C(sE — A) _1 B + D with a nonsingular matrix E. In this section, 
we will extend IRKA to descriptor systems. First, we establish the interpolatory H 2 
optimality conditions in the new setting. 

Theorem 4.1. Let G(s) = G sp (s) + P(s) be decomposed into the strictly proper 
and polynomial parts, and let G(s) = G sp (s) + P(s) have an r th -order strictly proper 
part G sp {s)j= C sp (sE sp - A sp ) _1 B sp . 

1. IfG(s) minimizes the H.i-error ||G — G||-h 2 over all reduced-order models with 
order strictly proper part, then P(s) = P(s) and G sp (s) minimizes the 



an r 



th 



%2-error ||G sp - G sp ||-H 2 - 
2. Suppose that the pencil sE sp — A sp has distinct eigenvalues {Ai}|_^. Let 
Yi and Zj denote the left and right eigenvectors associated with Xi so that 
A sp Zj =AiE sp Zi, y*A sp = Aiy*E sp , and y*E sp Zj = Then for q = C sp Zi 
and bj = y*B sp , we have 



G(-A,-)bi - G(-\i)b h cjG(-Xi) = cf G(-Ai), 
and cfG'(-Xi)bi = cjG'(-Xi)bi for i = 1, • • • , r. 



(4.3) 



Proof. 1. The polynomial part of G(s) and G(s) coincide, since, otherwise, the 
T^-norm of the error G(s) — G(s) would be unbounded. Then it readily follows that 
G sp (s) minimizes ||G sp (s) - G sp (s)||« 2 since G(s) - G(s) = G sp (s) - G sp (s). 

2. Since P(s) = P(s), the optimal model reduction problem for G(s) now 
reduces to the H.i optimal problem for the strictly proper transfer function G sp (s). 
Hence, the optimal H2 conditions of [TS] require that G sp (s) needs to be a bi-tan- 
gential Hermite interpolant to G sp (s) with {— Xi\ r i=1 being the interpolation points, 
and {ci} r i=l and {b^}[ =1 being the correspondin g left and right tangential directions, 



respectively. Thus, the interpolation conditions (4.3) hold since P(s) = P(s). □ 

Unfortunately, the H2 optimal interpolation points and associated tangent direc- 
tions are not known a priori, since they depend on the reduced-order model to be 
computed. To overcome this difficulty, an iterative algorithm IRKA was developed 
[T2l IT4] which is based on successive substitution. In IRKA, the interpolation points 
are corrected iteratively by the choosing mirror images of poles of the current reduced- 
order model as the next interpolation points. The tangential directions are corrected 
in a similar way; see [U [T3] for details. 

The situation in the case of descriptor systems is similar, where the optimal 
interpolation points and the corresponding tangential directions depend on the strictly 
proper part of the reduced-order model to be computed. Moreover, we need to make 
sure that the final reduced-model has the same polynomial part as the original one. 
Hence, we will modify IRKA to meet these challenges. In particular, we will correct 
not the poles and the tangential directions of the intermediate reduced-order model 
at the successive iteration step but that of the strictly proper part of the intermediate 
reduced-order model. As in the case of Theorem |3.1| the spectral projectors P; and 



P r will be used to construct the required interpolatory subspaces. A sketch of the 



resulting model reduction method is given in Algorithm 4.1 
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Algorithm 4.1. Interpolatory "H 2 optimal model reduction method 
for descriptor systems 

1) Make an initial selection of the interpolation points {o~i\ r i=1 and the 
tangential directions {bj}^ =1 and {ci}£ =1 . 

2) V/ = [ (triE - A)- 1 P / Bb 1 , . . . , (a r E - A)- 1 P / Bb r ] , 

W f = [ (o~%E - A)- T P^C T d, . . . , (<r r E - A)- T P^C T c r ] . 

3) while (not converged) 

a) A sp = WjAV /; Esp = WjEV /; B sp = W/B, and C sp = CV,. 

b) Compute A sp z; = AjE sp Zj and y*A sp = A 4 y*E sp with y*E sp Zj = 
where and Zj are left and right eigenvectors associated with . 

c) <Tj < \i, bf <- y*B sp and c { <- C sp Zj for i = 1, . . . ,r. 

d) V/ = [ (ciE - A) _1 P;Bbi, ((T r E — A) _1 P ; Bb r ] , 

W> = [ ( ffl B - A)- T P^C T c 1; . . . , (a r E - A)- T P^C T c r ] . 

end while 

4) Compute and smc/i i/iat ImfW^) = Im(I — P^) and 
Im(Voo)= Im(I-P r ). 

5) Set V = [V/, Voo] andW = [W/, W^. 

6) E = W T EV, A = W T AV ; B = W T B 7 C = CV, D = D. 



Note that until Step 4 of Algorithm |4.1[ the polynomial part isjiot included since 
the interpolation parameters result from the strictly proper part G sp (s). In a sense, 
Step 3 runs the optimal 'H.i iteration on G sp (s). Hence, at the end of Step 3, we 
construct an optimal %2 interpolant to G sp (s). However, in Step 5, we append the 
interpolatory subspaces with Voo and (which can be computed as described 

at the end of Section [3| so that the final reduced-order model in Step 6 has the 
same polynomial part as G(s), and, consequently, the final reduced-order model G(s) 
satisfies the optimality conditions of Theorem |4.1| One can see this from Step 3c: upon 
convergence, the interpolation points are the mirror images of the poles of G sp (s) and 
the tangential directions are the residue directions from G sp (s) as the optimality 
conditions require. Since Algorithm |4 . 1 1 uses the projected quantities P;B and CP,., 
theoretically iterating on a strictly proper dynamical system, the convergence behavior 
of this algorithm will follow the same pattern of IRKA which has been observed to 
converge rapidly in numerous numerical applications. 

Summarizing, we have shown so far how to reduce descriptor systems such that the 
transfer function of the reduced descriptor systems is a tangential interpolant to the 
original one and matches the polynomial part preventing unbounded Hco and H2 error 
norms. However, this model reduction approach involves the explicit computation 
of the spectral projectors or the corresponding deflating subspaces, which could be 
numerically infeasible for general large-scale problems. In the next two sections, we 
will show that for certain important classes of descriptor systems, the same can be 
achieved without explicitly forming the spectral projectors. 

5. Semi-explicit descriptor systems of index 1. We consider the following 
semi-explicit descriptor system 



E u xi(t) +E 12 x 2 (t) 


y(t) 



= AiiXi(t) + Ai 2 x 2 (t) + Biu(t), 
= A 2 ixi(t) + A 22 x 2 (t) + B 2 u(i), 
= Cixi(t) + C 2 x 2 (i)+Du(t) > 



(5.1) 
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where the state is x(i) = [xf(t), xJ(t)] T G R™ with xi(t) € R™ 1 , x 2 (t) € R" 2 and 
n x +n 2 = n, the input is u(t) G M m , the output is y(i) G M p , and E u , A n G M" lXni , 



2 



Ei 2 ,Ai2 G R niXn2 , A 2 i G IR" 2 *™ 1 , A 22 G ]R™ 2X " 2 , Bi G R niXm , B 
Ci G M pxni , C 2 G M px " 2 , D G R pxm . We assume that A 22 and En - EiaA^Aai 



are both nonsingular. In this case system (5.1) is of index 1. We now compute the 
polynomial part of this system. 



Proposition 5.1. Let G(s) be a transfer function of the descriptor system (5.1 1, 
where A 2 2 and En — E 12 A 22 1 A2i are both nonsingular. Then the polynomial part of 
G(s) is a constant matrix given by 



whe 



P(s) = C 1 M 1 B 2 + C 2 M 2 B 2 + D, 



Mi = (E u - E 12 A2 2 L A 21 )- 1 E 12 A 2 - 2 1 , (5.2) 
M 2 = — A 2 " 2 1 A2i(Eii - Ei 2 A 2 - 2 1 A 21 )- 1 Ei 2 A2 2 1 - a£. (5.3) 



Proof. Consider 
(sE - A) *B 
This leads to 



sEu - An sE i2 - A12 " 


-1 


Bi " 




' Fx( S ) " 


— A 2 i — A22 




B 2 . 




Fa(«) 



(sE u - An)Fi(s) + (sE 12 - A 12 )F 2 (s) - B b 
-A 21 F 1 (s)-A 22 F 2 (s)=B 2 . 



(5.4) 
(5.5) 



Solving ( 5.5 | for F 2 (s) gives F 2 (s) = -A 22 1 (B 2 + A 2 iFi(s)), and, thus, 

Fi(«) = ((aE u - An) - (SE12 - Ai 2 )A 22 1 A 2 i) _1 (Bi + (sE 12 - Ai 2 )A 22 1 B 2 ) 
implying that 

lim Fi(a) = (En - Ei 2 A 22 1 A 2 i) _1 E^A^Ba- 



Taking into account (|5.5|), we have 
lim P a (a) = 



- A 22 lA 21 ( E ll ~ Ei 2 A 22 1 A2l) E^A^ 1 - A22 1 



B 2 



Finally, note that P(s) = lim G(s) = lim (CiFjfs) + C 2 F 2 (s) + D), which leads to 

s— >oo s— ¥00 

the desired conclusion. □ 

We are now ready to state the interpolation result for the descriptor system 
(5.1). This result was briefly hinted at in the recent survey pQ. Here, we present it 



with a formal proof together with the formula developed for P(s) in Proposition 5.1 



As our main focus will be "T^-based model reduction, we will list the interpolation 
conditions only for the bi-tangential Hermite interpolation. Extension to the higher- 
order derivative interpolation is straightforward as shown in the earlier sections. 

Lemma 5.2. Let G(s) be a transfer function of the semi-explicit descriptor sys- 
tem (5.1). For given r distinct interpolation points {o~i} r i=1 , left tangential directions 
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<r and W 6 C nx ' be given 



{ci}i=i an d fight tangential directions {bi}[ =l7 let V £ 



V = [(tr 1 E-A) _1 Bbi J (ct>E - A) _1 Bb r ], 
W = [(cr x E- A)- T C T ci, (o>E- A)" T C T c r ]. 



(5.6) 
(5.7) 



Furthermore, let B and C be the matrices composed of the tangential directions as 

B = [ bi, . . . , b r ] and C = [ci, . . . , c r ]. (5.8) 
Define the reduced-order system matrices as 

E = W T EV, A = W T AV + C T D£, B = W T B C T f>, 

C = CV-DB, D = CiMiB 2 + C 2 M 2 B 2 + D. 



(5.9) 



Then the polynomial parts ofG(s) = C(sE — A) 1 B + D andG(s) match assuming^, 
is nonsingular, and G(s) satisfies the bi-tangential Hermite interpolation conditions 

G(* t )b t = Gfojbj, cj G(o-i) = cfG(a t ), cj G^b* = cf G'fo)^ 



/or i = 1, . . . , r, provided cr^E — A and 0"iE — A are both nonsingular. 



5.1 



we have 



Proof. Since E is nonsingular, lim G(s) = D. But by Lemma 

s— >oo 

D = lim G(s) ensuring that the polynomial parts of G(s) and G(s) coincide. The 

s— >-oo 

interpolation property is a result of [5J I22j , where it is shown that the appropriate 
shifting of the reduced-order quantities with a non-zero feedthrough term as done in 



(5.9) attains the original bi-tangential interpolation conditions hidden in V and W 



of (5.6 1 and (5.7 1, respectively. □ 



This result leads to Algorithm 5.1 which achieves bi-tangential Hermite inter- 



polation of the semi-explicit descriptor system (5.1) without explicitly forming the 
spectral projectors. 

Algorithm 5.1. Interpolatory model reduction for semi-explicit 
descriptor systems of index 1 

1) Make an initial selection of the interpolation points {<7j}£ =1 and the tangential 
directions {bj}£ =1 and {q}[ =1 . 

2) V= [(<7 X E- A)- 1 Bb 1 , (cr. r E- A)- 1 Bb r ], 

W= [(o-iE- A)- T C T c 1; (<j r E- A)- T C T c r ]. 



3) Define D = C1M1B2 + C 2 M 2 B 2 + D, where Mi and M 2 are defined in ( |5.2[ ) 
and (5.3), respectively. 



4) Define B = [ bi , . . . , b r ] and C = [ Ci , . . . , c r ] . 

5) E = W T EV ; A = W T AV + C T T)B, B = W T B - C T D, C = CV T>B. 



We want to emphasize that the assumption in Lemma 5.2 that E be nonsingular 
is not restrictive. This will be the case generically. If W and V are full-rank n x r 
matrices, the rank of the r x r matrix E = W T EV will be generically r as long as 
rank(E) > r. The fact that E will be full-rank is indeed the precise reason why we 



cannot simply apply Theorem 2.1 to descriptor systems 
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5.1. Optimal H2 model reduction for semi-explicit descriptor systems. 
Lemma |5.2| provides the theoretical basis for an IRK A-based iteration for "H 2 model 
reduction of semi-explicit descriptor systems. One naive approach would be the fol- 
lowing: Given system ( |5.1[ ), simply apply IRKA of [15 to obtain an intermediate 
reduced-order model G(s) = C(sE — A) _1 B + D. Of course, this will be generically 
an ODE and will not necessarily match the behavior of G(s) around s = 00. Thus, 
apply Lemma 5.2 to obtain the final reduced-model G(s) = C(sE — A) _1 B + D with 



E = E, A = A + C T D6, B = B C T D, 



C-DS, 



(5.10) 



where D is defined as in Lemma 5.2 While this shifting of the intermediate matrices 
by the D-term guarantees that the polynomial parts of G(s) and G(s) match, the "H 2 
optimality conditions will not be satisfied. The reason is as follows. Recall that the 
Ti.i optimality requires bi-tangential Hermite interpolation at the mirror images of the 
reduced-order poles. The intermediate model G(s) satisfies this but since it does not 
match the p olyno mial part, the resulting % 2 error is unbounded. Then constructing 
G(s) as in (5.10), we enforce the matching of the polynomial part but G(s) still 
interpolates G(s) at the same interpolation po ints a s G(s), i.e., at the mirror images 
of the poles of G(s). However, clearly due to (5.101, the poles of G(s) and G(s) are 
different; thus G(s) will no longer satisfy the optimal %2 necessary conditions. In 
order to achieve both the mirror-image interpolation conditions and the polynomial 
part matching, the D term modification must be included throughout the iteration, 



not just at the end. This results in Algorithm 5.2 



Algorithm 5.2. IRKA for semi-explicit descriptor systems of index 1 

1) Make an initial shift selection {<7i}; =1 and initial tangential directions {bi}f =1 
and {cj}J" =1 . 

2) V r = [(01E- A)- x Bbi, (a- P E- AJ^Bbr], 
W r = [ (<7tE - A)- T C T c l7 . . . , (o>E - A)- T C T c r ] . 



3) Define D = C1M1B2 + C 2 M 2 B 2 + D. where Mi and M 2 are defined in ( |5.2| 
and (5.3), respectively. 

4) Define B = [ bi, . . . , b r ] and C = [ ci , . . . , c r ] . 

5) while (not converged) 

a) E = W T EV,_A = W T AV + C T T)B, B = W T B - C T D, C = CV T>B. 

b) Compute Y*AZ = diag(Ai, . . . , A r ) and Y*EZ = I r , where the columns of 
Z = [zi,...,z r ] and Y = [yi,...,yy] are, respectively, the right and left 
eigenvectors of AE — A. 

c) <7j < \i, bj 4— y*B and q 4— Czi for i = 1, . . . , r. 

d) V= [(aiE- A)- 1 Bbi, (cr r E - A)- x Bb r ] , 

W = [ (<riE - A)- T C T c l7 . . . , (a r E - A)- T C T c r ] . 
end while 

6) E = W T EV, A = W T AV + C T T)B, B = W T B - C T T), C = CV - T>B. 



The next result is a restatement of the above discussion. 
Co roll ary 5.3. Let G(s) be a transfer function of the semi- explicit des criptor 

A B 



5.2 



The 



system (5.1) and let G(s) = C(sE — A) B + D be obtained by Algorithm 
G(s) satisfies the first-order necessary conditions of the "H 2 optimal model reduction 
problem. 
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5.2. Supersonic inlet flow example. Consider the Eulcr equations modelling 
the unsteady flow through a supersonic diffuser as described in [3T]. Linearization 
around a steady-state solution and spatial discretization using a finite volume method 
leads to a semi-explicit descriptor system ( 5.1 1 of dimension n = 11730. For simplicity, 



we focus on the single-input single-output subsystem dynamics corresponding to the 
input as the bleed actuation mass flow and the output as the average Mach number. 

It is important to emphasize that applying balanced truncation to this model 
is far from trivial because of difficulty of solving the Lyapunov equations. Instead, 



we apply the proposed method in Algorithm 5.2 to obtain an "H 2 -optimal reduced- 
model of order r = 11, where the only cost are sparse linear solves and the need for 
computing the spectral projectors are removed. As pointed out in [21], the frequencies 
of practical interest are the low frequency components. Figure [5T2] shows the amplitude 
and phase plots of G(iu>) and G(iu>) for uj E [0, 25] illustrating a very accurate match 
of the original model. The resulting model reduction errors are 
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Fig. 5.1. Supersonic inlet flow model: amplitude and phase Bode plots of G(s) and G(s). 



6. Stokes-type descriptor systems of index 2. In this section, we consider 
a Stokes-type descriptor system of the form 

Enii(t) = A u xi(t) + Ai2X 2 (t)+B 1 u(t) ) 

= A 2 ixi(i)+B 2 u(i), (6.1) 
y(t) - Cixi(t) + C 2 x 2 (i) + Du(t), 

where the state is x(t) = [xf(t), *2(t)] T E E" with x x (£) E R ni , x 2 (t) E M" 2 and 
n 1 +n 2 = n, the input is u(t) E R m , the output is y(i) E R p , and En, An E R niXn \ 
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Ai 2 € R niXn2 , A 2 i € M" 2Xni , Bx e R" lXm , B 2 G R" 2Xm , Ci G M px ™\ C 2 G M px ™ 2 , 
and D G R pxm . We assume that E n is nonsingular, A 12 and A^ h ave both full 
column rank and A2iE^ 1 1 Ai2 is nonsingular. In this case, system (6.1 ) is of index 2. 

In [T7|, the authors showed how to apply ADI-based balanced truncation to sys- 
tems of the form (6.1) without explicit projector computation. Here, we extend this 



analysis to interpolatory model reduction and show how to reduce ( |6.1[ ) optimally 
in the % 2 - n orm without computing the deflating subspaces. Unlike |17j . En is not 
assumed to be symmetric and positive definite, and A21 is not assumed to be equal 
to Af 2 . 



First, consider system (6.1 1 with B 2 = 0, as the case of B 2 =^ follows similarly. 



Following the exposition of [17] , consider the projectors 

II, = I - E u 1 Ai 2 (A 2 iE 11 1 Ai 2 )- 1 A 2 i, 
n r = I- A 12 (A 21 E n 1 A 12 )- 1 A 21 E 1 - 1 1 . 



Then the descriptor system (6.1) can be decoupled into a system 



IIjEuIIr*i(t) = lIiAnllrXiCtJ+njBiuCO 
y{t) = Cn r xi(t) + Du(t) 



(6.2) 



with 



C = Ci - C 2 (A 2 iE u 1 Ai2) _1 A 2 iE n 1 Aii, 
D = D C 2 (A 21 Er 1 1 A 12 )- 1 A 21 Er 1 1 B 1 , 

and an algebraic equation 

x 2 (t) = -(A2iE n 1 A 12 )- 1 A 2 iE n 1 A 11 x 1 (i) - (A2iE n 1 A 12 )- 1 A 21 E u 1 B 1 u(t). 

By decomposing II; and II r as 



with e^-, & r . 3 G r«i x ("i-«2) suc h that 



and defining xi(t) = r2 xi(f), system (6.2 1 becomes 



©f 2 E u e rl ii(t) = 0f 2 A n r)lXl (i) + 0f 2 B lU (t) 
y(t) = C0 ril xi(f) + Du(i), 



(6.3) 



(6.4) 



(6.5) 



Then the reduction of the descriptor system (6.1) is equivalent to the reduction of 



system (6.2) or (6.5). However, the beauty of this equivalence lies in the observation 
that the matrix ;2 En0 rl is nonsingular. Therefore, standard model reduction 



procedures for ODEs can be applied to system (6.5), and the obtained reduced-order 



model will approximate the descriptor system (6.1 ). It is important to emphasize that 



even though (6.2) and (6.5) are equivalent to (6.1), the ultimate goal of this section is 



to develop an interpolatory model reduction method that does not require the explicit 
computation of either the projectors 11;, II r or the basis matrices 0z, 2 , r ,i- For this 
purpose, define the matrices 



£ = n,E u n r , A = n ; Ann r , s = n ; Bi, e = cn r . (6.6) 
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In interpolation setting, the matrix of interest will be er£ — A with a G C. Luckily, 
several key properties of £ + tA, t G C, were already introduced in [17] . However, 
we present these results in terms of <r£ — A instead of £ + tA. 

Lemma 6.1. Let ®i2 and ® r ^i be the matrices defined in (6.3) and let a G C be 
such that a®J 2 E'ii® r 1 — 0f 2 An0 r 1 is nonsingular. The matrix defined as 

(<r£ - A) 1 := e Pil (ffe£ 2 E u e ril - e^AuQ^)- 1 ©^ (6.7) 

satisfies 

(cr£ - AY(a8. - A) = n r and (cr£ - A)(aE - A) 1 = II,. 
Similarly, the matrix defined as 

(a£ T - A 7 ) 1 := ©^(a®^,®^ - ©^A^e^)" 1 ©^ (6.8) 

satisfies 

(cr£ T ~A T ) I (a8. T -A T ) =IlJ and (<r£ T - A T )(a£ T - A T )' = . 



Proof. Following a similar argument to that in |17j , the proof of the first equality 
follows directly from (6.3) and (6.7). Indeed, we have 

(aE-AY(aE - A) = ril (a6[ 2 E„0 rjl - 0^A 11 r a)" 1 ©Ln/( ( TE 11 - A u )n r 

= e ril ((7e£ a E u e ril - ©J^®^)- 1 ®^^ - A n )® rA ®r, 2 

= ®r,l® T r,2 = 

The remaining equalities follow similarly. □ 



At first glance, the definition of the generalized inverses in (6.7) and (6.8) may 
seem to be irrelevant for model reduction of the descriptor system (6.1 ). Recall that 
reducing (6.1) is equivalent to reducing system (6.2) and the interpolatory projection 
method for (6.2) will require inverting (<t£ — A) and (<r£ T — A T ). However, these 
inverses do not exist. As a result, definitions (|6.7[ ) and (6.8) become pivotal in order 
to achieve interpolatory model reduction of (6.2) and, thereby, of (6.1) as shown in 
the next theorem. 

Theorem 6.2. Let s — er, // G C be such that the matrices 

s0^ 2 Eu0 r)1 - 0[ 2 Au0 ril and sW T E n V - W T A n V 

are invertible. Define the reduced- order model 

G(s) = CV(sW T E n V - W T AnV)- 1 W T Bi + D. (6.9) 
Let b G C m and c G C p be fixed nontrivial vectors. 

1. //(ffE-ill'Sbe Im(V) C lm(0 r ,i) and (/i£ T - A T ) I e r c G Im(W) C lm(0i, 2 ), 
then G(cr)b = G(a)b and c T G(n) = c T G(fi). 

2. //, in addition, a = /i, then c T G'(<r)b = c T G'(a)b. 

Remark 6.1. Before presenting the proof, we want to emphasize that this interpo- 



lation result is different than the usual interpolation framework given in Theorem 2.1 
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where the projection matrices V and W are constructed using A, E, B and C and 
then the projection is applied to the same quantities. In Theorem |6.2[ however, the 



projection matrices V and W are constructed using the system matrices of (6.2 1 , 
namely A,&,15 and C. But then the projection (model reduction) is applied to the 
system matrices of (6.1), namely En,An,Bi and C. Thus, the proof will serve to 
fill in this important gap. 



Proof. Since systems (6.1) and (6.5) are equivalent, they have the same transfer- 
function given by 



G(s) = C& r As&T 2 B n @ ril - ©^Au©^)- 1 ©^! + D. 



Since ©; 2 Eii© r 1 in (6.5) is nonsingular, we make use of Theorem 
W such that 



2.1 



Define V and 



V = © r .iV 



and 



W = 0/ 2 W. 



(6.10) 



Pluging these matrices into (6.9), we obtain that 

G(s) = C0 r ,iV(sW T 0i r 2 Eii0,. a V - W T 0^ 2 A 1 i© I . i iV)- 1 W T 0^ r 2 B 1 + D. 



Moreover, it follows from (6.4) that V = 0^ 2 V and W = ©^W. To prove the first 
claim in part 1, we note that (6.3) implies that 



(6.11) 



Since (cr£ - 7l) J S b e Im(V ), there exists qel r such that (<r£ - A/Sb = Vq. 
Using (6.7) (6.10) and ( |6.11 1, this equation can be written as 

©^(ae^En©^ - e^AnB^j-'e^Bib = © r ,iVq. 

The left multiplication by 0^ 2 gives 

K 2 En0 r! i - e^AuQ^O-^^Bib = Vq. 

Hence, (o-Q^En©^ - &l 2 A 11 & rl )~ 1 &f 2 B 1 b E Im(V). Then it follows from 
that G(er)b = G(<r)b. The equation c T G(a) = c T G(cr) ned 



2.1 



Theorem 

similarly. The proof of part 2 follows from part 3 of Theorem |2.1| □ 

It should be noted that the conditions Im(V) C Im(® r .i) and Im(W) C lm(0; 2 ) 
in part 1 of Theorem |6.2| are automatically fulfilled if for given interpolation points 
{ <J i}i=i) {A*i}i=i an d tangential directions {bj}[ =1 , {q}[ =1 , we choose 

Im(V) = span{(cri£ - A) 7 2bi, . . . , (cr r £ - A) z Sb r }, 
Im(W) = span{(^i£ T - A T ) / C T Ci, . . . , (a r £ T - A T ) I & r c r }. 

6.1. Computational issues related to the reduction of index-2 descrip- 
tor systems. Even though Theorem 6.2 shows how to enforce interpolation for the 
descriptor system ( |6.1[ ), the spectral projectors are still implicitly hidden in the defi- 
nitions of (cr£ — A) 1 and (cr£ T — A T ) ! . It has been shown in [T7] how to compute 
the matrix-vector product (£ + rA) I f for a given vector f without explicitly forming 
(£ + tA) 1 . This approach can also be used in interpolatory model reduction, where 
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the quantities of interest are (c£ — A) I 'Bb and (/i£ T — A T ) I G T c. The proof of the 
following result is analogous to those in [T7], and, therefore, it is omitted. 

Lemma 6.3. Let s = a,fi be such that s9f 2 Eu0 r 1 — ®J 2 A-n® r ± is invertible. 
Then the vector 

v=(aE-A) I 'Bb (6.12) 

solves 



o-En - An A12 " 




V 




' Bib " 


Aai 




z 








(6.13) 



and the vector 

w = (fiE T -A T ) I e T c (6.14) 

solves 



V&Ti - A-u A^ " 




w 






Af 2 




q _ 








(6.15) 



From a computational perspective of implementing Theorem 6.2 the importance 



of this result is clear. To achieve interpolation, Theorem 6.2 relies on computing the 
quantities (cr£ — A) 1 lib and (er£ T — A T ) I G T c 1 both of which involve the computation 



of 0/ 2 and r i. However, Lemma 6.3 illustrates that the computation of these basis 



matrices is unnecessary and only the linear systems ( 6.13 1 and (6.151 need to be solved 



This observation leads to Algorithm 6.1 below for interpolatory model reduction of 
Stokes-type descriptor systems of index 2. 



Algorithm 6.1. Interpolatory model reduction for Stokes-type 
descriptor systems of index 2 

1) Make an initial selection of the interpolation points {(Ti}|* =1 and the tangent 
directions {bi}[ =1 and {Ci\l =1 . 

2) For i = 1, . . . , r, solve 



OiEn — An 


A12 ' 








" Bib, " 


A 2 i 







z 

















" C T c 4 " 


A T 

A 12 







q 








3) V = [vi, . . .,v r ],__ 

4) E = W T EnV, A 



W [wi,.. 
= W T AnV, 



B = W T B 1; C = CV, D = D. 



Once a computationally effective bi-tangential Hermite interpolation framework is 
established for index-2 descriptor systems, extending it to optimal Ti.2 model reduction 
via IRK A is straightforward and given in Algorithm 6.2 It follows from the structure 
of this algorithm that upon convergence the reduced model G(s) = C(sE— A) _1 B+D 
satisfies the first-order conditions for Ti.2 optimality. 

Remark 6.2. As shown in [17j . the general case B2 7^ can be handled 
similar to the case B2 = 0. First note that the state ~x.i(t) can be decomposed as 
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Algorithm 6.2. IRK A for Stokes-type descriptor system of index 2 

1) Make an initial shift selection {cr,}^ =1 and initial tangent directions {b;}^ =1 and 
{C|! 

2) Apply Algorithm 



6.1 



to obtain E, A, B, C and D. 

3) w/iiZe (not converged) 

a) Compute Y*AZ = diag{\\ 1 . . . , A r ) and Y*EZ = I, where the columns of 
Z = [zi,...,z r ] and Y = [yi,...,y r ] are, respectively, the right and left 
eigenvectors of AE — A. 

b) <7j < Ai, b^" -f- y*B and Cj <— Cz^ /or i = 1, . . . , r. 



c) Apply Algorithm 
end while 



6.1 



to obtain E, A, B, C and D. 



xi(i) = xo(i) + x 9 (t), w/iere x g (t) = -E 11 1 Ai 2 (A 2 iE 11 1 A i2 ) 1 B 2 u(i) and x (i) sa- 
tisfies A 2 iXo(t) = 0. After some algebraic manipulations, this leads to 



njEiin r io(t) = n z Ann r xo(t) + njBu(i), 

y(t) = Cn r x„(t)+Du(^C 2 (A 21 E u 1 A 12 )- 1 B 2 u(i), 



(6.16) 



where 



C = d - C 2 (A 2 iE n 1 Ai 2 ) _1 A 2 iE n 1 Aii, (6.17) 
B = Bi — A 11 E n 1 A 12 (A 21 E u 1 A 12 )- 1 B 2 , (6.18) 
D = D C 2 (A 21 E n 1 A 12 )- 1 A 21 E n 1 B 1 . (6.19) 

Therefore, the B 2 ^ case extends to the interpolation framework as well by defining 

£ = n ; Eiiii r , A = n;Aiin r , 23 = n ; B, e = cn r 



and applying Theorem 6.2 with £, A., 15, C and 2) = D — ,sC 2 (A 2 iE 11 1 Ai 2 ) X B 2 
instead o/£, A, 25, C and 2). 

6.2. Numerical results for Oseen equations. The model borrowed from [17] 
is obtained by discretizing the Oseen equations and describe the flow of a viscous and 
incompressible fluid in a domain n e R 2 representing a channel with a backward 
facing step. A spatial dis creti zation using the finite element method leads to the 
index-2 descriptor system O) with En, An <= m 5520x5520 , A^A^ € R 5520 * 761 , 



Bi G R 552uxb ; g 2 g K 7bix£r-^ i g ^2x552^ ^ £ R Jx7bi ) D = , see |T7| for more 
details on the model. Note that B 2 7^ and the transfer function grows unbounded 
around s = 00. 

We approximate this system by a model of order r = 20 using the balanced trun- 
cation method as described in [17) and the T~Li optimal model reduction method given 



in Algorithm 6.2 The amplitude Bode plots of the full model and two reduced-order 



models depicted in Figure 6T clearly illustrate that interpolation-based Algorithm |6.2| 



leads to a high-fidelity reduced model replicating the full-order transfer function with 
almost no loss of accuracy and matching the performance of the balanced trunca- 
tion method. The accuracy of this interpolation-based method is due to the fact 
that we do not choose the interpolation points in an ad hoc fashion; instead Algo- 
rithm |6.2| iteratively leads to 'H 2 optimal interpolation points. As the difficulty in 
computing "H 2 norm of the error is clear, we approximately compute the relative "Hoo- 
error ^ G « P G ^p^ 00 for both reduced-order models by sampling the imaginary axis. 
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These errors for the balanced truncation method and Algorithm [6^2] are, respectively, 
3.3284 x 10~ 6 and 8.9663 x 10 -6 . Both reduced-order models are highly accurate. It is 
expected that the 'H 00 -error in balanced truncation will be smaller than that in IRKA. 
While our method tries to minimize the T^-norm, the balanced truncation method is 
tailored towards reducing the "Hoo-norm. Indeed, these numbers are further signs for 
the success of the interpolatory-based model reduction method as it produces a very 
accurate model, almost matching the accuracy of the balanced truncation approach. 
These observations are similar to those on IRKA whose Hoo-norm behavior was close 
to or even better in some cases than that of balanced truncation [TJ [TS] . 



Amplitude Bode Plots 



Fig. 6.1. Oseen equation: amplitude Bode plots of the full and reduced models 
To further illustrate the accuracy in the reduced-order model computed by Al- 



gorithm 6.2 we display the time domain response plots resulting from two different 
input selections. In the left pane of Figure |6.2| we plot the outputs for the input 
selections Uj(t) = sin(6ii) for i = 1, . . . , 6 (recall that the system has 6 inputs). The 
figure illustrate a perfect match between the outputs of the full and reduced-order 
systems. Error in the outputs for the same input selection is given in the right pane 



of Figure 6.2 Note the difference in the scale of the error plot compared to the actual 
output; the error is four orders of magnitude smaller. We repeat the same experi- 
ments with Uj(i) = sin(zi) for i = 1, ... ,6 and reach the same conclusions as shown 
in Figure [6~3] 



7. Conclusions. For interpolatory model reduction of descriptor systems, we 
have introduced subspace conditions that not only guarantee interpolation conditions 
but also automatically enforce matching the polynomial part of the transfer function, 
thus preventing the error grow unbounded. We have also extended the optimal H2 
interpolation point selection strategy to descriptor systems. For the index- 1 and 
index-2 descriptor systems, we have shown how to construct the reduced-order models 
without computing the deflating subspaces corresponding to the finite and infinite 
eigenvalues explicitly. Several numerical examples have supported the theoretical 
discussion. 
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Output 1 foru.(t) = sin(6 i t) 



Error in Output 1 for u.(t) = sin(6 i t) 
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Fig. 6.2. Oseen equation: (left) time domain response for Ui(t) = sm(6it); (right) error 
in time domain response for Ui(t) — sin(6it). 



Output 1 for u.(t) = sin(i t) 



Error in Output 1 for u.(t) = sin(i t) 
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Fig. 6.3. Oseen equation: (left) time domain response for Ui(t) — sin(it); (right) error 
in time domain response for Ui(t) = sin(ii). 
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